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Abstract We investigate the secular dynamics of two-planet coplanar systems 
evolving under mutual gravitational interactions and dissipative forces. We con- 
sider two mechanisms responsible for the planetary migration: star-planet (or 
planet-satellite) tidal interactions and interactions of a planet with a gaseous disc. 
We show that each migration mechanism is characterized by a specific law of 
orbital angular momentum exchange. Calculating stationary solutions of the con- 
servative secular problem and taking into account the orbital angular momentum 
leakage, we trace the evolutionary routes followed by the planet pairs during the 
migration process. This procedure allows us to recover the dynamical history of 
two-planet systems and constrain parameters of the involved physical processes. 

Keywords Secular evolution • Migration ■ Dissipation laws ■ Tidal force ■ Drag 
force 

1 Introduction 

Recently, Michtchenko and Rodriguez (2011) have purposed a new method for a 
qualitative study of the planet migration originated by a generic dissipative mech- 
anism. It has been shown that, under assumption that the dissipation processes are 
sufficiently slow, the evolutionary routes of migrating planets in the phase space 
are traced by stationary solutions of the conservative secular problem. Therefore, 
the modeling of planet migration consists in the calculation of families composed 
by Mode I and Mode II stationary solutions, parameterized by the mass ratio, for 
all possible values of planetary semi-major axes and orbital angular momentum of 
the system. 
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It is presently accepted that the main mechanisms responsible for planet migra- 
tion are the star-planet (and planet-satellite) tidal interactions and gravitational 
interactions of a planet with a gaseous disc. Dissipative interactions remove (or 
add) orbital energy from the planetary system, producing changes of semi-major 
axes. In the case of tidal interactions, orbital energy is transformed into thermal 
energy, which is dissipated in the interior of the tidally deformed body. Migration 
of the system can be either toward or outward the central body, depending on the 
properties of the tidal interactions and some parameters such as Love numbers and 
quality factors (Darwin 1880; Jeffreys 1961; Goldreich and Soter 1966; Hut 1981; 
Dobbs-Dixon et al. 2004; Ferraz-Mello et al. 2008). During disc-planet interactions, 
orbital energy is exchanged with a surrounding gaseous protoplanetary disc and, 
generally, the planetary orbital decay occurs (for a review, see Armitage 2010 and 
references therein). In addition, migration can also occur through gravitational 
scattering between planets and a remnant planetesimal debris (e.g. Fernandez and 
Ip 1984; Kirsh 2009). 

Dissipative forces also alter the orbital angular momentum of the system. In- 
deed, tidal torques modify the rotational angular momenta of the interacting bod- 
ies, which are transferred to the orbits of planets or satellites (e.g. Mignard 1979; 
Correia et al. 2008). During disc-planet interactions, the angular momentum is 
exchanged through gravitational torques between the planets and the disc (e.g. 
Lin et al. 1996; Goldreich and Sari 2003; Masset et al. 2006). In Michtchenko and 
Rodriguez (2011), we have purposed a generic approach to describe the exchange 
of the orbital angular momentum of the system during migration. We have intro- 
duced the orbital angular momentum leakage, defined as a portion a of the orbital 
angular momentum variation produced by migration (i.e. expansion or shrinkage 
of the planetary orbits), which is extracted (or added to) from the planet system. 

In this paper, we show that each physical process is characterized by a specific 
law of orbital angular momentum leakage. We show that a is a function of the 
eccentricity of the planet affected by the dissipative force. Moreover, we show 
that a also depends on a set of physical parameters related to planets, star and 
disc. Finally, owning the characteristic a-function, we construct the evolutionary 
routes of the system for each specific mechanism driving the planet migration. This 
approach provides us with a general idea of how the system could evolve under a 
variety of migration conditions and physical models. 

In Section 2 we briefly introduce some basic aspects of the secular dynamics 
of two-planet systems evolving under dissipative forces. Section 3 addresses tidal 
interactions, where we include the cases of a close-in planet and a satellite interact- 
ing with their parent star and planet, respectively. Migration driven by disc-planet 
interaction is discussed in Section 4. The dissipation is modeled through a drag 
Stokes-like force. Finally, Section 5 is devoted to conclusions. 



2 Some aspects of the secular dynamics under dissipation 

2.1 Orbital angular momentum exchange 

We consider a three-body system consisting of a central star with mass mp and 
two coplanar planets with masses mi and m2- We assume that the system is far 
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away from any mean-motion resonance. Hereafter, the indices i = 1,2 stand for 
the inner and outer planets, respectively. 

In the astrocentric reference frame, the orbital angular momentum of the sys- 
tem is given, up to second order in masses, as 

Lorb = m'ly/aiil - el) + m2y/a2{l - e|), (1) 

where = rrn^J G{mo + rrii), ai and are planetary semi- major axes and eccen- 
tricities, and G is the gravitational constant. 

The migration of planets is originated by dissipative forces which affect the 
energy and the orbital angular momentum of the planets and results in variations 
of their semi-major axes and eccentricities. If dissipation is sufficiently slow (its 
rate is smaller than the proper frequency of the secular motion), the long-term 
variations of the orbital elements can be separated into two components: one is 
produced by secular interactions with the companion and the other is due to 
external interactions. Hence, we can write 

• sec I • ext /r\\ 

di = a-i + a,i , (2) 

p. _ p. sec ,. ext 

where the indices "sec" and "ext" stand for secular and external contributions of 
the total variation of each element. 

The secular theory provides that '''"^ = (to first order in masses) and, as a 
consequence, ai = di We stress that, since we are studying the secular evolution 
of the system, short-period variations of the orbital elements (of the order of orbital 
periods) are omitted in the analysis. Thus, Equations (Hl-Q should be considered 
as the averaged equations describing the long-term variations of semi-major axes 
and eccentricities. In addition, the contribution of resonant terms is also neglected. 

Assuming that the planetary masses are unaltered during migration, the time 
variations of Lorbi in terms of the time variation of the planetary semi-major axes, 
ai, and eccentricities, e^, is written as 

. m[y/l- ej , mi^/oTei . m'a^/l - e| . 7712^0262. 
Lorh = 5-7= cii ei + -— as €2. (4) 

The total angular momentum of the whole system (including the three-body 
system under study and the external component) is conserved during migration. 
Thus, we can write that Loj-b + ^ext = const, where the external component of 
the total angular momentum is denoted as Lext- It should be emphasized that the 
generic definition "ext" does not necessarily mean that the exchange of the orbital 
angular momentum occurs with an exterior medium. For instance, in the case of 
tidal interactions of planets with a central star, the 'external' contribution to the 
angular momentum comes from the rotation of the star, as described in the next 
section. On the other hand, in the case of disc-planet interactions, there is a flux 
of angular momentum between the planet system and an external protoplanetary 
disc, as described in Section [H 

For sake of simplicity, we restrict our investigation to the case in which only one 
planet is directly affected by dissipative forces and, consequently, its semi-major 
axis is changed. This implies that the semi-major axis of the other planet remains 
unchanged during secular evolution of the system. Assuming that the i body 
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is affected by dissipative forces and combining Equations Q - (|4}, we use the 
conservation of the total angular momentum to find that 




(5) 



and 



. ext_ (ijlfi). , ViEfif (f^\ 

where i ^ j and ^ 0. Equation ([5]) shows that, under mutual secular perturba- 
tions, the eccentricities of the planetary orbits oscillate in anti-phase. Equation ([5} 
also shows that the evolution of ej depends implicitly on the migrating evolution 
of the companion, through the external variations of and e^. 
Equation ([6} can be re-written as 

e,^'^* = (l-a)%l^a„ e, ^ 0, (7) 



a^ ^ 0, (8) 



where 

2(1-1 -^ext 
a = : — , 

with Lj being the partial angular momentum of the body affected by dissipation. 
The function a can thus be used as a convenient measure of the non-conservation 
of Lorb, since it quantifies the variation of the external component of the orbital 
angular momentum of the system. Indeed, if a = 0, Lext remains constant and L^rh 
is conserved during migration. According to Equation dTJ, the limit case = 
implies that ctj = for a = 0, that is, the migration must be stopped when the 
planet orbit is circularized in the system with no exchange of angular momentum. 

For a 7^ 0, the portion 1 — a of the I/orb- change produced by hi is absorbed by 
the system through the variation of ej, while the rest (a) is transferred to Lext. 
We say that there is a leakage (or gain) of the orbital angular momentum in the 
system. According to Equation ([7}, for a = 1, the dissipative force produces no 
change of the eccentricity and the change in Lorb due to hi is totally transferred 
to Lext. This particular case of the angular momentum exchange was used to 
modeling of planet-planetesimal interactions (Malhotra 1995). 

The range of theoretical values of a is large; according to Equation ([8}, a 
can be positive or negative, depending on the signs of hi and Lext- Moreover, ct 
is a function of the orbital elements of the migrating planet and, consequently, 
varies during migration. Knowing that each physical process is characterized by 
a specific law of angular momentum exchange, in the next sections, we develop 
the Q-function for several kinds of dissipative interactions. Using Equation ([7]), we 
write a in general form as 

Q = 1 - -^(em par) \ , (9) 

where J^(ej,par) is a characteristic function of the migration process, defined 
through the condition 

^=.F(e„par)^, (10) 

Gi (Xi 

where the vector par is composed of physical parameters of the process under 
study. 
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2.2 Evolutionary routes in the phase space 

As shown in Michtchenko and Rodriguez (2011), under assumption that dissipation 
is weak and slow, the evolutionary routes of the migrating non resonant planets are 
traced by the Mode I and Mode II stationary solutions of the conservative secular 
problem (see also Hadjidemetriou and Voyatzis 2010). The ultimate convergence 
and the evolution of the system along one of these secular modes of motion are 
determined by the condition that the dissipation rate is smaller than the proper 
secular frequency of the system. We have shown that, for values of a less than 
1, the Mode I secular solution (characterized by aligned orbits) plays a role of 
an attractive center when planetary orbits diverge and the condition ^ a\la.2 < 
m2/m\ is satisfied, or, when planetary orbits converge and ^Joijaa > m2/mi. In 
opposite cases, the Mode II (characterized by anti-aligned orbits) is the attractive 
center of the migrating two-planet system. 

In practice, the calculation of evolutionary tracks is simple. The current loca- 
tion of the system in the phase space provides the starting values of the orbital 
angular momentum, iorb, stud the semi-major axis of the planet directly affected 
by dissipation. The set, composed by Lorbj masses and semi-major axes, uniquely 
defines two stationary solutions (e^^, ej) of the conservative secular problem at the 
current configuration. In order to model migration process, the semi-major axes 
of the affected planet is incremented by Aai (which can be either positive or nega- 
tive), while the semi-major axis of the other planet is kept unaltered. The orbital 
angular momentum of the system is then corrected by the amount 

Ml 

where the function a is defined by Equation ([9} . For a new set of and -Lorb values 
(with fixed masses), we calculate new values of stationary solutions, (61,62). The 
procedure is repeated until the system reaches domains with no possible stationary 
solutions. The obtained values of (6^,62) are then plotted on the representative 
planes (ni/n2, e^), where rii (i = 1, 2) are mean motions of the planets. The family 
of stationary solutions shows two evolutionary routes, one of which the system will 
follow during migration. 

In the described above process, stationary solutions can be calculated using the 
precise semi-analytical approach developed in Michtchenko and Malhotra (2004), 
with no restrictions on the values of planet eccentricities. An alternative is the use 
of the expansion of the disturbing function in series of axjan-, ex and 62, given by 

N I N' / \ ' 

= X! ^i,mA,i'^ ( ^ 1 eSz' cos(mZiti7), 

where _R; „j ;/ are numerical coefficients which do not depend on the physical and 
orbital parameters of the system and Am = tui — taJ2 is the difference of planetary 
longitudes of pericenter. 

The main difficulty in the calculation of migration routes is that the charac- 
teristic function J^(ei,par), which is needed to obtain a through Equation Q, is 
unknown a priory. Thus, the next sections of this paper will be devoted to calculate 
the function T for some specific dissipative processes. 
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3 Tidal interactions 

We first consider tiie case of tidal interactions, which affect orbital elements and 
rotations of the deformed bodies, while energy is dissipated due to internal friction. 
The tidal force is inversely proportional to the 7 power of the distance between 
interacting bodies and is effective in both planetary and satellite systems (e.g. 
Mignard 1979). In this paper, we assume that only the central and close-in bodies 
are mutually affected by the raised tides, whereas the outer companion is not 
directly affected by the tidal force. The orbital planes of planets are assumed to 
be coincident. 



3.1 A system with a close-in planet 

We consider a short-period planet orbiting a slow-rotating central star (Qq ni), 
where Qq and ni are the angular velocity of rotation of the star and the mean 
orbital motion of the planet, respectively. 

The long-term variations (averaged over the inner planet orbital period) of the 
planetary elements, due to the combined effects of tides raised on the star and on 
the inner planet, are written, up to ©(ef), as: 

< ai > = -|nia^^s[(l + 23e?) + 7efZ)] (11) 



and 



where 



with 



< ei""* >= ^|nieiaj;'^s[9-h7D], (12) 



D = p/2s, (13) 



9 fci mo „5 J . 9 fco mi 5 , , 

p=-- 7?i and s=-- i?o (14 

2 Qi mi 4 Qo "lo 

being the strengths of stellar and planetary tides, respectively (Dobbs-Dixon et 
al. 2004; Ferraz-Mello et al. 2008, Rodriguez and Ferraz-Mello 2010). fe, and Qi 
are the Love number and the dissipation function of the deformed body, where 
the indices i = and 1 stand for the central and inner bodies, respectively. In 
above equations, it is assumed a quasi-synchronous rotation state of the close-in 
planet {Qi ~ ni). In addition, is is also assumed a linear dependence between 
phase lags and the corresponding frequencies appearing in the decomposition of 
the tidal potential of each body (Mignard 1979; Ferraz-Mello et al. 2008). For sake 
of simplicity, we consider that Qi are constant; however, it should be kept in mind 
that their dependence on frequencies can be important in the long-term evolution. 
The reader is referred to Efroimsky and Williams (2009) for more discussions on 
the frequency dependence of Qi. 

Combining Equations pop and (|lip - p2p . we obtain the characteristic function 

9 -h 7-D 

^^^^'^^ = 2[l + (23 + 7i3)ef] ' (^'^ 
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Fig. 1 The Lo^b ^leakage a as a function of ei , in the case of a close-in planet tidal evolution. 
The curves are parameterized by several constant values of D. 



where par is defined by the parameter D. The Loj.)-,— leakage or the function a is 
then calculated through Equation Q to give 



1 



(9 + 7D)ef 



(l-ef)[l + (23 + 7D)e?] 



(16) 



Figure [T] shows the variation of a with ei. We observe that, for all values of D, 
a varies in the range between and 1. The immediate consequence is that the tidal 
decay of a close-in planet is generally accompanied by the loss of orbital angular 
momentum of the system. The portion 1 — q is absorbed by the system, damping 
the eccentricity of the inner planet. The amount a is transferred to the central 
star and accelerates the stellar rotation (see discussion of the next paragraph) . 

The stellar rotation changes due to tides on the star and, from the above 
discussion, it is clear that its variation plays a role of external source of the angular 
momentum. The rate of variation of the rotational angular momentum is given by 
iext = irot = Cof2o, whcrc Co is the stellar moment of inertia (assuming that the 
inner planet has reached stationary rotation, we can neglect the spin variation of 
the planet; e.g. Rodriguez et al. 2011 for more details). Hence, using the definition 
of a given in Equation ([8]), we obtain, for di ^ 0, 



2Coai 
Liai 



On 



(17) 



Since ai < and a > during orbital decay of the inner planet, above equation 
results in SIq > 0, as previously discussed. 

On one hand, Equation (|17p shows that, for a = 0, is a constant and there is 
no tides on the star (i.e. s = 0). In this case, the orbital angular momentum of the 
system is conserved. Inversely, knowing that D = p/2s, s = means that D — oo 

= 0.5e^^. Hence, to second order in ei 
0. Note that, for a ~ (very large D), 



and, according to Equation p5|) . T{ei,D) 
and according to Equation ([9]) , we have a = 
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Fig. 2 Evolutionary routes of CoRoT-7 system parameterized by constant values of D. The 
angular momentum exchange is determined by the a-function ||9jl defined through II16II . The 
current position of the system is shown by a star symbol, whereas the locations of some mean- 
motion resonances are indicated by dashed lines. A schematic view of the migration scenario 
is shown in the box illustration at the top panel. 

the migration is halted when the orbit circularizes, that is, ai becomes constant 
because there is no orbital angular momentum transfer and ei = 0. 

On the other hand, when planetary tides are neglected, we have p = D = 
and, up to second order in ei, (1 — a) = 9ei/(l + 23ef ). This last expression shows 
that, after circularization of the inner planet orbit (ei = 0), the orbital angular 
momentum change due to the orbital decay is totally transferred to the stellar 
spin (a = 1). Evidence for excess of rotational angular momentum of the parent 
stars due to the tidal interaction with close-in giant planets have been already 
addressed in previous studies (e.g., Pont 2009). 

It is important to note that, during tidal interactions, a is not constant but 
varies as a function of ei, according to Equation (|16p . The sign of a is always 
positive, meaning that the decreasing of the orbital angular momentum due to the 
orbital decay of the inner planet is compensated by both the damping of ei and 
the increase of Hq. 

Figure[2]shows the evolutionary tracks of a migrating pair of planets. They were 
obtained as discussed in Section 12.21 where the orbital angular momentum varia- 
tion for each curve was defined by the function a, according to Equation (|16p . An 
application to the CoRoT-7 planetary system, characterized by very small current 
eccentricities (Ferraz-Mello et al. 2011) was considered, whose results are plotted 
in the planes {n\/n2, e,j). The adopted migration scenario is schematically shown 
on the top of the figure. The divergent migration, together with the condition 
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raijmx > ^Ja\/a2, results in the Mode I {A-oj = 0) of secular motion as an attrac- 
tive center (Michtchenko and Rodriguez 2011). The current position of the system 
is marked with a star symbol and the arrows indicate the direction of the evolution 
(orbital decay of the inner planet). 

We can note that migration tracks are sensible to the values of the parameter D, 
and consequently, to the orbital angular momentum exchange. For small D (a ~ 1), 
corresponding to the nearly-total loss of the angular momentum variation produced 
by the orbital decay of the inner planet, both eccentricities vary only slightly 
during evolution. In this case, when the strength of the stellar tide is dominating, 
the orbital configuration of the system in the past would be not different from the 
present one. 

The situation is very different in the case of large values of D, when the strength 
planetary tide is dominating. Large D-values imply that the orbital angular mo- 
mentum of the system is almost conserved (a ~ 0). In this case, the orbital config- 
uration of the system in the past should be characterized by high eccentricities, as 
higher as larger are D-values. Therefore, in the cases with dominating planetary 
tides, the origin of such high eccentric planetary orbits must be investigated. 

Note that, during migration, the two-planet system could cross several mean- 
motion resonances, resulting in the temporary excitation of eccentricities but with- 
out trapping, since the orbits are divergent. Several numerical simulations per- 
formed in Rodriguez et al. (2011) and Michtchenko and Rodriguez (2011), show 
that the system ultimately returns to the stationary secular solution after leaving 
the mean-motion resonance. 



3.2 A system of satellites 

We consider now the case of a pair of satellites orbiting a rapidly rotating parent 
planet {fio » ni). We assume that only the inner satellite is affected by tidal 
interactions with the planet. According to the linear tidal model, the averaged 
variations of the orbital elements are given by 

< ai > = |nia7*s[(l + (27/2 - 7D)e?] (18) 

and 

< er* >= |nieia^5s(ll/2-7D) (19) 
o 

(Goldreich and Soter 1966; Ferraz-Mello et al. 2008). Note that the only difference 
with respect to the previous case is in the tides raised on the planet (or central 
body). We still note that the mean variations can be positive or negative, depend- 
ing on the value of D. Hence, the combined effect of planet and satellite tides can 
induce either inward or outward migration, as well as either damping or excitation 
of eccentricity. The sign of * depends on the value of D in such a way that 
^ext > Q if < 11/14. The sign of a\ depends on both D and ei-values. For small 
eccentricities (ei < 0.2) and moderate D (in the range from 1 to 5), we would ex- 
pect outward migration and damping of eccentricity, a typical behavior observed 
in the tidal evolution of Solar System satellites (e.g. Peale 1999). In this case (i.e., 
(ii > and e^"* < 0), the total angular momentum conservation implies that a > 1 
(see Equation ([7])) and, according to Equation p7p . flo < 0. The Earth-Moon tidal 
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Fig. 3 Dependence of a on ei , in logarithmic scale, for the tidal evolution of a close-in satellite. 
Illustrations arc shown for different values of D. Here, a can adopt values larger than unity, 
but not in the close-in tidal evolution case (see Figure.^ for comparison). 



evolution, with the satellite moving away from the Earth and the increasing length 
of day of the planet, corresponds to this scenario. 



Combining Equations ©-([TOl) and (IT8l)-(IT9l), we have 



and 



_1 {11/2 -7D)el 

{l-e?)[l + (27/2-7D)e?]- ^ ^ 

Figure |3] shows the variation of a with ei, for different values of D. With the 
exception of the example with D = 0, a > 1 always, indicating that the external 
source of angular momentum (planet rotation) injects angular momentum in the 
satellite system, which is used to expand the inner planet orbit and circularize 
both inner and outer orbits. 



Figure. H] shows the migration routes of a hypothetical system in which the 
central body is an Earth-like planet, the inner satellite has the Moon mass and 
the outer companion is three times smaller than the Moon. The curves are param- 
eterized by several values of D. For D in the range from to 5, the inner satellite 
migrates outward when the eccentricity of its orbit is not too large (ei < 0.2, 
according to Equation (HH}). Since the semi-major axis of the outer satellite is 
unaltered during the secular interaction with its companion, we have a convergent 
migration and, together with the condition m2/mi < a\/a2, results in the Mode 
I {Azu = 0) as the attractive center. The direction of migration along evolutionary 
routes is shown by arrows in Figure. H) 
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Fig. 4 The same as in Figure [2] except the angular momentum exchange is determined by 
the a— function defined through I I21I I. 

For D > 1, the satellite eccentricities decrease during migration: larger is the 
value of D, more rapid is the damping of eccentricities. In this way, the satellite 
system approaches to main mean-motion resonances with low-eccentricity orbits 
which allow a smooth capture into one of these resonances (e.g. Tittemore and 
Wisdom 1988). The behavior of the system described by D = is different, be- 
cause in this case the eccentricities of the satellites increase during migration (see 
Equation (|19p ). As a consequence, the capture inside a low-order mean-motion 
resonance would be strongly improbable. 

4 Disk-planet interactions 

In this section we study the migration of a two-planet coplanar system assum- 
ing that dissipative forces affect only the outer planet. We also assume that the 
two planets are interacting secularly, that is, the planets are far enough from 
any mean-motion resonance. Disc-planet interactions are the natural mechanism 
which results in the orbital decay of the outer planet due to energy and angular 
momentum exchange with an outer gaseous disc (e.g. Kley 2000, 2003). 

4.1 Evolutionary tracks under Stokes interactions 

For sake of simplicity, we suppose that the outer planet is forced to migrate under 
the action of a dissipative drag force (Stokes-like force) of the type 
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f = -10 "(^2 - 7V2c), 



(22) 



where v > Q and V2c is the Keplerian circular velocity at the astrocentric distance 
r-2. Due to the negative radial pressure gradient, the gas velocity is a bit less than 
the circular velocity at the same point. Consequently, the value of the factor 7 
should be slightly smaller than unity (see Adachi et al. 1976; Patterson 1987). 
7 is also frequently used as a free parameter of the problem, in order to model 
different conditions of disc-planet interactions (e.g. Beauge et al. 2006). Several 
previous works have investigated the evolution of a two-planet system evolving 
under a Stokes drag force (e.g., Beauge at al. 2006; Hadjidemetriou and Voyatzis 
2010), with the majority devoted to the study of the evolution inside mean-motion 
resonances. 

In order to obtain the variations of the semi-major axis and the eccentricity 
of the outer planet produced by the force f, we use the Euler-Gauss's equations 
(Brouwer and Clemence 1961): 



da2 
IF 

de2 
~dt 

dzU2 

dt 

dq 

dt 



n2 



[fr 62 Sinu2 + /t(l + 62 COSU2)] 



(23) 



n2 a2 



fr s'mU2 + — f 1 + 62 COS U2 

62 v 



n2 0,2 62 



-/r COSU2 + /t 



r2 



02(1 



r2 

0.1 

-h 1 I sin U2 



62 



[(02(1 - el) COSU2 - 2e2r2)/r - (r2 + 02(1 - 6^)) sin 1^2 /t] 



where 1% is related to the mean anomaly, ^2, through I2 = '2 +/ 'f^idt. fr and /( are 
the radial and transverse components of f , whereas U2 is the true anomaly of the 
outer planet. The first-order averaging over the outer planet orbital period gives, 
to fourth order in 62 and assuming m2 <^ mo: 



< a2> = -10" 

< 62 > = -10" 
<ZU2> = 



''a2[2(l - 7) + 7(56i/8 + 1196^/512)] 



"^7 62(1 - 1361/32) 



(24) 
(25) 
(26) 



(for comparison, see Beauge and Ferraz-Mello 1993). It follows from above equa- 
tions that the disc-planet interactions produce orbital decay and circularization. 
If we set 7=1 (see Hadjidemetriou and Voyatzis 2010) and 62 = into Equation 
(|24p . then we have a2 = 0, indicating that there is no orbital decay of the outer 
planet in the case of circular orbit. 

Using Equations ©-([TOl) and (I2l|)-(l25l), we obtain 



-^(62,7) 



7(l-13ei/32) 



2(1 - 7) -h 7(562/8 + 11964/512) 



and 



1 - 



27(l-13el/32)el 



(1 - el) [2(1 - 7) + 7(56^/8 + 1196|/512)] 



(27) 
(28) 



where par is given by the parameter 7. 
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Fig. 5 Variation of a with 62, for two values of 7, corresponding to the case of disc-planet 
interaction simulated with a drag force of the type |(22j. Note that, depending on 62, o. can 
adopt negative values (see text for details). 



Figure [5] shows the function 0(62,7) parameterized by two values of 7. We note 
that, as in the case of star-planet tidal interactions (see Sect. 13. ip . q is always 
smaller than 1, which means that the orbital decay of the outer planet is always 
accompanied by damping of the eccentricity of its orbit (see Equation (O). We 
also see that, for e2 = 0, q = 1 for any value of 7 (we have discussed that 7 < 1), 
indicating that, after circularization of the planet orbit, the change of Lorb due to 
the orbital decay is totally transferred to Lext, in this case, to the gaseous disc. 

In contrast with the case of star-planet tidal interactions, the function a is 
not saturated at (see Figure [T|) , but decreases monotonously with eccentricities 
towards negative values. According to Equation ^ , negative values of a imply in 
the loss of angular momentum by the disc during the orbital decay of the outer 
planet (Lext < 0). Thus, during migration, the pair of planets can either gain or 
lose orbital angular momentum, depending on the eccentricity of the migrating 
planet. 

A transition between leakage or gain of angular momentum occurs at a = 0, 
when 62 reaches a critical value e'^ which is determined from Equation (|28p , for a 
given 7. For instance, for 7 = 0.995, e'^ ~ 0.085. For 62 = e", the orbital angular 
momentum of the planet system is conserved (since a = 0), however, this condition 
is only transitory. For ei < e" , we have < q < 1, indicating that a portion a 
of the orbital angular momentum change due to migration is removed from the 
system and transferred to the disc (Lext > 0, see Equation (H])), while the portion 
1 — a is absorbed by the system altering the value of e2. When 62 > e" , we have 
a < and the system gains angular momentum from the disc. Note that, in the 
latter case, the planet system gains angular momentum, despite the semi-major 
axis of the outer planet is decreasing. 

The evolutionary routes of the planet pair, obtained for several values of the 
parameter 7 and for the angular momentum exchange defined by (|28p , are shown 
in Figure (6] The hypothetical system is composed of a Sun-like central body and 
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two planets of equal mass. The migration scenario shown on the top of the figure 
corresponds to the outer planet moving towards the inner planet. In this case, 
migration is convergent and the Mode II plays a role of an attractive center, since 
m2/mi > axjai (see Michtchenko and Rodriguez 2011). Planet eccentricities 
are damped during migration (see Equation (|25p l. while the system approaches to 
one of the main mean-motion resonances, whose locations are indicated in Figure 
|6] Because of small eccentricities, the probability of capture in a mean-motion 
resonance would be high, mainly for larger 7-values. It is worth noting that, to 
achieve the current position of the system (marked by a star symbol), a past orbital 
configuration with higher eccentricities is needed for all considered values of 7. 

It is important to stress that, during the migration path, the system crosses 
several mean-motion resonances, as shown in Figure [6] One should keep in mind 
that, when the rate of dissipation is small, the system may be trapped one of such 
resonances. Since we restrict our investigations to the secular evolution of the sys- 
tems, we fix the initial configurations of planet pairs far away from mean-motion 
resonances and, according to our model, neglect the possibility of a resonant cap- 
ture. 
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4.2 The factor K 



Several works have simulated the disc-planet interaction adopting a simplified 
model of planetary migration (e.g. Lee and Peale 2002; Beauge et al. 2006, among 
others). The method avoids the use of hydrodynamic numerical simulations in- 
volving an accurate description of the interaction process between the gas and the 
planet body. The model introduces the relationship between the migration rates 
of semi-major axis and eccentricity as follows 



(29) 

where if is a constant (positive) parameter. Given a value of K, the migrating 
systems are modeled through numerical integrations of the equations of planetary 
motion, with constant perturbations in the orbital elements according to Equation 
(1291). 

As shown in Beauge and Ferraz-Mello (1993) and Gomes (1995), the time 
evolution of orbital elements due to the action of the drag force given by (|22p is, 
to first order in eccentricity: 



a2{t) = a2oexp{—At) and e2(i) = 620 exp(— _Et), 

where A and E are the inverse of the e- folding times of 02 and 62, respectively. In- 
troducing above equations into Equation (|29p we obtain that K = E/A. Moreover, 
using Equations (|24p - ([25)) we also obtain that, to first order in 62, K = 7/2(1 —7). 
In this way, we have shown that the condition (|29|1 . with a constant value of K, can 
be directly obtained from the averaged variations of the orbital elements produced 
by the Stokes force. It is worth noting that, for 7 = 0.995, which is the usually 
adopted value (Beauge and Ferraz-Mello 1993), we obtain K = 100, which matches 
the value determined empirically by Lee and Peale (2002) for the GJ 876 resonant 
system (the planets b and c evolve inside the 2/1 mean-motion resonance). 

The relationship (|29p is frequently applied in moderate and even high eccentric- 
ity domains, specially when the evolution of resonant pairs of planets is modeled. 
However, comparing Equations (|10p and (|29p we note that K = T{e2,'y), that is, 
the assumption of constant K is not appropriate in the case of eccentric orbits and 
the function J^(e2,7), given by Equation (|27p . should be used instead K = const 
(see also Kley et al. 2004). The variation of if as a function of 62 is shown in 
Figure [3 



4.3 A general drag force 

The Stokes force (|22p belongs to the class of a more general dissipative force given 
by 

f = -10-''p(r2)(v2-7V2c), (30) 

where p{r2) is the density profile of the gas in the disc (see Smart 1960). In this 
work, we assume that p{r2) = , with /? > (see Kominami and Ida 2002, for 
an example with /3 = 2). Some works also adopt a more general dependence on 
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0.1 0.2 0.3 
62 



Fig. 7 Variation of K with 62, for two values of 7. When £2=0 and 7 = 0.995, A" = 100 is 
obtained, reproducing the result of previous works (Lee and Peale 2002). Moreover, K strongly 
varies with £2 and for high eccentric orbits is almost independent on the value of 7. 



the relative velocity, introducing powers of (v2 — 7V2c) larger than 1. This case is 
related to high Reynolds numbers and correspond to regions of turbulence in the 
gas medium. 

Using the Euler-Gauss equations with the external force (|3Up and applying the 
averaged procedure (over orbital periods), we obtain the long-term variations of 
semi-major axis and eccentricity of the outer planet, up to 0{e\), as follows 

<a2 >=-10-'^4-^[2(l-7) + eigi(7,/3)+e|g2(7,/?)], (31) 
<e2 >=-10-'^a-^ 62 [7(1-/3) + /3 + el e;3(7,/3)], (32) 

where 

2^1 = 3/? + + 7 (1 - 2/3 - , 
1652 = 7/3 -h y/32 + 5/33 i/3^ 

403 = -3/3 + |/3^ + i/33 + 7(/? _ 1) _ _ . 

It is clear that the case /3 = corresponds to the previously studied Stokes force 
(see Equations (f24)) - ([251) 1. 

Interesting features can be highlighted from the above equations. On one hand, 
we note that fixing 02 and 62, the rate of migration increases when /J < 1 and 
decreases if /? > 1. Moreover, for /3 = 1, < 02 > does not depends on 02. On the 
other hand, the rate of 62 - damping decreases for /3 > and it is independent of 
02 for the Stokes drag (/3 = 0). 
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Fig. 8 Variation of a with 62 for the case of a general dissipative force given by II30II . Three 
values of the power /3 are illustrated for 7 = 0.995. 



Through Equations (|9ll- (fT0)l and (f3T |l - (f32|l . the characteristic function J^(e2, 7, /?) 
is given by 

T{. ^ R\- [7(l-/?)+/? + eig3(7,/?)] .oox 
>le2,7,W [2(l-7) + eiei(7,/3) + e|e2(7,/3)] ^ ^ 

and the function a 

„ ^ 1 _ 2ei[7(l-/3) + /3 + eig3(7,/3)] 

(l-e2)[2(l-7) + eiei(7,/J) + 4e2(7,/?)]' 

where par is given by two parameters of the system, 7 and /3. Figure [8] shows 
the variation of 0(62,/?, 7) as a function of £2, for 7 = 0.995. We show curves 
parameterized by /3 = 0,1,2. Note that a has a strong dependence on 62 for all 
/3. For high eccentricity (e2 > 0.2), a < for all values of /3 illustrated. When 
a < 0, the planetary system gains orbital angular momentum, because a2 < and, 
looking at Equation (|8]), we have Lext < 0. Moreover, the gain is larger as smaller 
is the power fi. In the limit of very small eccentricity (ei < 0.05), a is almost 
independent on /3. 

The dependence of the planet evolution on the power /3 is illustrated in Figure|9l 
where we show the migration paths of the fictitious two-planet system, previously 
analyzed in Section. 14. II Several values of /? are used, with 7 fixed at 0.995. We can 
note that, at least for the given 7, the different density profile distributions produce 
only quantitative differences in the damping of eccentricities, which are significant 
only at high eccentricity domains. The convergent migration of the planets can 
results in the crossing or capture in a mean-motion resonance, as in some previous 
discussed cases. 
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A = 0, A < 

^ -o *-© 




5 e 7 8 9 




n / n 2 

Fig. 9 Evolutionary routes for a disc-planet interaction simulated through the general force 
OOll . parameterized by constant values of the power /3 and 7 = 0.995. The angular momentum 
exchange is determined by the corresponding a— function, given by Equation I I 3411 

5 Conclusions 

In this work we model analytically the orbital angular momentum exchange of two- 
planet coplanar systems evolving under dissipative forces. Two migration mech- 
anisms were considered: tidal interactions in star-planet and planet-satellite sys- 
tems and gaseous disc-planet interactions. Our approach is based on the model 
developed in Michtchenko and Rodriguez (2011), which shows that the angular 
momentum exchange between the orbital and the exterior components of the total 
angular momentum can be calculated through the a-function, referred to as the 
orbital angular momentum leakage. For each dissipative mechanism considered, a 
was calculated as a function of the planet eccentricity of the migrating planet and 
physical parameters involved in the process. 

Using the obtained a-function, stationary solutions of the secular conservative 
problem were calculated for each specific migration scenario and for several val- 
ues of physical parameters. Stationary solutions provide evolutionary routes that 
the system would follow in the process of planetary migration and allows us to 
understand the dynamical history of the system evolving under dissipative forces. 

The tidal interactions between a close-in planet with its host star results in the 
orbital decay of the inner planet. The angular momentum exchange between the 
orbital component and the stellar rotation imposes constraints on the q- function in 
such a way that < q < 1. During migration, the planet system always loses some 
part of the orbital angular momentum, which is used to accelerate the rotation of 
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the star. Large (small) values of the parameter D (ratio between the strengths of 
planetary and stellar tides) are associated to weak (strong) stellar tides, enabling 
a substantial conservation (dissipation) of the orbital angular momentum of the 
system. 

In the planet-satellite tidal interaction, the situation is generally opposite. In 
this case, the migration of the inner satellite is predominantly outward. For typical 
values of the parameter D (in the range between 1 and 5) and moderate eccentric- 
ities (ei < 0.2), the a-function is always larger than 1. This means that angular 
momentum is injected into the satellite system, with the planet spin as the supplier 
source. 

Disk-planet interactions modeled through a drag Stokes-like force, have shown 
that Of < 1 always. Angular momentum can be removed (0 < a < 1) or injected 
(a < 0) into the planet system. We have shown that the factor K (the ratio 
between eccentricity damping and orbital decay), which is frequently adopted as 
a free parameter of the dissipative problem, is a function of the outer planet 
eccentricity. In addition, K also depends on the parameter 7, related to the disc 
properties. Moreover, the assumption of constant K is only valid in the domain of 
very small orbital eccentricities. 

Finally, the development of the a-function for each dissipative process, and 
the consequent calculation of evolutionary routes allow us to reassemble the start- 
ing configurations and migration history of the planet systems on the basis of 
their current orbital configurations. In addition, the analysis of the orbital angu- 
lar momentum evolution during migration of the system allows us to constrain 
parameters of the involved dissipative physical processes. 
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